
load('grid_100.mat')
N_temp = 10000;

d2s_dzdz = NaN(N_temp,1);
d2s_dzdx = NaN(N_temp,1);
ratio_temp = NaN(N_temp,1);

for i=1:N_temp

    [~,temp] = max(data.z(i,:));

    temp_diff = setdiff(1:J,temp);

    ind_good2 = temp_diff(1);

    argument = normcdf([data.x(i,temp) data.x(i,temp_diff) data.z(i,temp) data.z(i,temp_diff)]);

    d2s_dzdz(i) = fun_d2sigma_new(theta_NPD(:,ss),argument,m_regr,combos_all_unique,J+1,J+2)*normpdf(data.z(i,temp))*normpdf(data.z(i,ind_good2));

    d2s_dzdx(i) = fun_d2sigma_new(theta_NPD(:,ss),argument,m_regr,combos_all_unique,J+1,2)*normpdf(data.z(i,temp))*normpdf(data.x(i,ind_good2));

    ratio_temp(i) = d2s_dzdz(i)/d2s_dzdx(i);
    
end

[num_temp_ord,ind_num] = sort(d2s_dzdz);
[den_temp_ord,ind_den] = sort(d2s_dzdx);
ratio_of_trimmed_means_second(ss) = mean(num_temp_ord(N_temp*0.25:N_temp*0.75))/mean(den_temp_ord(N_temp*0.25:N_temp*0.75));

[ratio_temp_ord,ind] = sort(ratio_temp);
trimmed_mean_of_ratios(ss) = mean(ratio_temp_ord(N_temp*0.25:N_temp*0.75));


% first derivatives

ds_dz = NaN(N_temp,1);
ds_dx = NaN(N_temp,1);
ratio_temp = NaN(N_temp,1);

for i=1:N_temp

    argument = normcdf([data.x(i,:) data.z(i,:)]);

    ds_dz(i) = fun_dsigma(theta_NPD(:,ss),argument,m_regr,combos_all_unique,J+1)*normpdf(data.z(i,1));

    ds_dx(i) = fun_dsigma(theta_NPD(:,ss),argument,m_regr,combos_all_unique,1)*normpdf(data.x(i,1));

    ratio_temp(i) = ds_dz(i)/ds_dx(i);
    %i
end

[num_temp_ord,ind_num] = sort(ds_dz);
[den_temp_ord,ind_den] = sort(ds_dx);
ratio_of_trimmed_means_first(ss) = mean(num_temp_ord(N_temp*0.25:N_temp*0.75))/mean(den_temp_ord(N_temp*0.25:N_temp*0.75));

display('*** RATIO OF TRIMMED MEANS -- SECOND ***')
ratio_of_trimmed_means_second(ss)

display('*** RATIO OF TRIMMED MEANS -- FIRST ***')
ratio_of_trimmed_means_first(ss)

display('*** TRIMMED MEAN OF RATIOS -- SECOND ***')
trimmed_mean_of_ratios(ss)


